# Goal: Create Figures 1 & A5
# by Jennifer Pan and Yiqing Xu

## Correlations of Six Dimensions

rm(list=ls(all=TRUE))

question.list <- list(paste0("s1_",1:14), # poli
                      paste0("s2_",1:14), # econ
                      paste0("s3_",1:14), # nati
                      paste0("s4_",1:7), # trad
                      paste0("s5_",1:7), # equi
                      paste0("s6_",1:7)) # ethn

qs <- unlist(question.list)


##########################################
## Matrix plot: Sample 1 (Figure 1a)
##########################################

library(psych)
library(corrplot)
library(GGally)
col <- RColorBrewer::brewer.pal(9, "PuBu")[8]

d <- haven::read_dta("data/sample1.dta")
names(d)

###### Index ######
M <- as.data.frame(d[,paste0("index_",c("poli","econ","nati","trad","equi","ethn"))])
for (i in 1:ncol(M)) {M[,i] <- as.numeric(M[,i])}

colnames(M) <- c("Politically\nLiberal","Pro-Market","Nationalistic", 
                 "Traditionalist","Pro-Social\nEquality","Accommodating\nMinorities")

## Raw distribution
colnames(M) <- c("Politically Liberal","Pro-Market","Nationalistic",  
                 "Traditionalist","Pro-Social Equality","Accommodating Minorities")
pdf("graphs/corr_index_q.pdf", height = 10, width = 10)
p <- ggpairs(M,lower = list(continuous = wrap("points", alpha = 1, size=0, colour = col)),
             upper = list(continuous = wrap("cor", size = 5, colour = 1)))  +
  theme(axis.text = element_text(size = 14, vjust = 1, color = "black"),
        axis.title =element_text(size=18, face="bold")) + theme_bw()
print(p, left = 0.2, bottom = 0.2)
graphics.off()


##########################################
## Matrix plot: Sample 2 (Figure 1b)
##########################################


d <- haven::read_dta("data/sample2.dta")
names(d)

###### Index ######
M <- as.data.frame(d[,paste0("index_",c("poli","econ","nati","trad","equi","ethn"))])
for (i in 1:ncol(M)) {M[,i] <- as.numeric(M[,i])}

## Raw distribution
colnames(M) <- c("Politically Liberal","Pro-Market","Nationalistic",  
  "Traditionalist","Pro-Social Equality","Accommodating Minorities")
pdf("graphs/corr_index_s.pdf", height = 10, width = 10)
p <- ggpairs(M,lower = list(continuous = wrap("points", alpha = 1, size=0, colour = col)),
             upper = list(continuous = wrap("cor", size = 5, colour = 1)))  +
    theme(axis.text = element_text(size = 14, vjust = 1, color = "black"),
          axis.title =element_text(size=18, face="bold")) + theme_bw()
print(p, left = 0.2, bottom = 0.2)
graphics.off()



##########################################
## Scree Plot (Figure A5)
##########################################

#library(lavaan)
#library(Amelia)
library(psych)

## Sample 1
d <- as.data.frame(haven::read_dta("data/sample1.dta"))
names(d)
d1 <- d[, qs]
dim(d1)
d1$id <- 1:nrow(d1)
set.seed(123)
out <- scree(d1[,1:63])
fv <- out$fv
pcv <- out$pcv


pdf("graphs/scree_q.pdf")
par(mar = c(5, 4, 1, 1))
plot(1, type = "n", ylim = c(-0.5, 7), xlim = c(1,63),
     xlab = "", ylab = "", axes = FALSE)
mtext("i'th Principal Component or Factor", 1, 2.5, cex = 1.5)
mtext("Eigenvalues of Components and Factors", 2, 2.5, cex = 1.5)
abline(h = 1, col = "#AAAAAA50", lwd = 4, lty = 1)
abline(h = seq(-1,20,0.5), col = "#AAAAAA50")
abline(v = seq(0,64,by=2), col = "#AAAAAA50")
lines(fv, col = "gray30", lwd = 2, lty = 3)
points(fv, pch = 16, col = "white", cex = 1.8)
lines(pcv, col = 1, lwd = 2)
points(pcv, pch = 16, col ="white", cex = 1.8) 
points(fv, col = "gray30", pch = 1, cex = 1.1)
points(pcv, pch = 16, col =1, cex = 1.2)
axis(1, at = seq(0,60,10), cex.axis = 1.3)
axis(2, at = 0:10, cex.axis = 1.3)
legend("topright", c("Principal Component", "Factor"), pch = c(16, 1),
       col = c(1, "gray30"), bty = "n", cex = 1.5)
box()
graphics.off()


## Sample 2
d <- as.data.frame(haven::read_dta("data/sample2.dta"))
names(d)
d1 <- d[, qs]
dim(d1)
d1$id <- 1:nrow(d1)
set.seed(123)
out <- scree(d1[,1:63])
fv <- out$fv
pcv <- out$pcv

pdf("graphs/scree_s.pdf")
par(mar = c(5, 4, 1, 1))
plot(1, type = "n", ylim = c(-0.5, 7), xlim = c(1,63),
     xlab = "", ylab = "", axes = FALSE)
mtext("i'th Principal Component or Factor", 1, 2.5, cex = 1.5)
mtext("Eigenvalues of Components and Factors", 2, 2.5, cex = 1.5)
abline(h = 1, col = "#AAAAAA50", lwd = 4, lty = 1)
abline(h = seq(-1,20,0.5), col = "#AAAAAA50")
abline(v = seq(0,64,by=2), col = "#AAAAAA50")
lines(fv, col = "gray30", lwd = 2, lty = 3)
points(fv, pch = 16, col = "white", cex = 1.8)
lines(pcv, col = 1, lwd = 2)
points(pcv, pch = 16, col ="white", cex = 1.8) 
points(fv, col = "gray30", pch = 1, cex = 1.1)
points(pcv, pch = 16, col =1, cex = 1.2)
axis(1, at = seq(0,60,10), cex.axis = 1.3)
axis(2, at = 0:10, cex.axis = 1.3)
legend("topright", c("Principal Component", "Factor"), pch = c(16, 1),
       col = c(1, "gray30"), bty = "n", cex = 1.5)
box()
graphics.off()



